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Abstract 

In this work, single-asperity contact mechanics is investigated for positive and negative work of adhesion Ay. In the latter case, 
finite-range repulsion acts in addition to hard-wall constraints. This constitutes a continuum model for a contact immersed in a 
strongly wetting fluid, which can only be squeezed out in the center of the contact through a sufficiently large normal load F$. As 
for positive work of adhesion, two stable solutions can coexist in a finite range of normal loads. The competing solutions can be 
readily interpreted as contacts with either a load-bearing or a squeezed-out fluid. The possibility for coexistence and the subsequent 
discontinuous wetting and squeeze-out instabilities depend not only on the Tabor coefficient |ij but also on the functional form of 
the finite-range repulsion. For example, coexistence and discontinuous wetting or squeeze-out do not occur when the repulsion 
decreases exponentially with distance. For positive work of adhesion, the normal displacement mainly depends on F^, Ay, and |ij 
but - unlike the contact area - barely on the functional form of the finite-range attraction. The results can benefit the interpretation 
of atomic force microscopy in liquid environments and the modeling of multi-asperity contacts. 



Introduction 

The continuum description of single-asperity contact mechanics 
has received much attention in the last few decades. This is in 
large parts because it describes force-displacement curves rather 
accurately down to nanometer scales relevant to atomic force 



microscopy (AFM) [1-3]. The contributions to the linear elas- 
ticity of (frictionless) single-asperity contacts most central to 
this work are the following: Hertz [4] solved the contact 
mechanics of a parabolic tip pressed against a substrate for 
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hard-wall repulsion. He found that the contact area A c and the 
separation between the two solids d both disappear continuous- 
ly with F^ 3 as the normal load squeezing the two solids 
together, F N , approaches zero. Derjaguin, Muller, and Toporov 
(DMT) [5] included adhesion as a long-range interaction, in 
which case adhesion effectively acts as a normal load that is 
independent of the contact-radius a c . Johnson, Kendall, and 
Roberts (JKR) [6] solved the problem for short-range adhesion, 
for which the attractive surface forces act directly at the contact 
line. Unlike Hertz and DMT, JKR found finite values for A c and 
d at pull off. Tabor [7] introduced a dimensionless parameter, 
now known as Tabor coefficient, allowing one to estimate if 
a contact is closer to DMT or to JKR theory. He actually recog- 
nized that DMT and JKR describe the opposite limits of long- 
and short-range forces, respectively. This had not been known 
before but was soon confirmed in numerical simulations by 
Muller, Yushenko, and Derjaguin [8]. Lastly, Maugis [9] used a 
cohesive-zone model introduced by Dugdale (MD) and found 
analytical solutions for intermediate-range adhesion at arbitrary 
values of 

Although single-asperity, linearly-elastic, adhesive contacts 
mechanics is a rather mature field [10], two key issues remain 
worth addressing: First, only few studies have considered the 
case of negative work of adhesion [1 1,12], Ay < 0, specifically 
finite-range repulsion between two surfaces acting in addition to 
hard-wall repulsion. In particular, the concept of the Tabor 
coefficient has not yet been applied to negative Ay. Therefore, I 
investigate if there are different regimes for Ay < 0 as is the case 
for contacts with Ay > 0, which are classified as JKR for large 
jix and as DMT for small This includes a study of which 
parameters determine the behavior near squeeze-out as well as a 
comparison to the behavior near pull-off for Ay > 0. For the 
latter, it is straightforward to deduce from established results 
how the « c (^n) relation depends on the Tabor coefficient in the 
DMT and the JKR limit. Specifically, a c - a p oc (F N + F p ) K for 
- ~Fp> where F p and a p are pull-off force and pull-off 
radius, respectively. They both depend on jij just like the expo- 
nent k, e.g., fl p > 0 and k = 1/2 in the JKR limit, while a p = 0 
and k = 1/3 for DMT. In the present comparison of squeeze-out 
(finite-range repulsion) versus pull-off (finite-range attraction), 
I also study whether the exponent k changes continuously 
between JKR and DMT or discontinuously - as assumed 
implicitly in the Carpick-Ogletree-Salmeron (COS) model [1]. 

The second motivation for this paper is that it has not yet been 
investigated sufficiently how the (precise) functional form for 
adhesive interactions affects contact mechanics - assuming that 
all continuum parameters, from normal load to Tabor coeffi- 
cient, are identical. It is only established that there is little sensi- 
tivity in the limits of large and zero Tabor coefficients. Yet, 



when studying contact-mechanics between macroscopic, adhe- 
sive, rough surfaces in the context of continuum-mechanics, one 
would want to know how to best reach the JKR limit, which 
appears to be the relevant limit for that application. For 
example, it is used implicitly in Persson theory for nominally 
flat, adhesive surfaces [13]. In fact, the current work was initi- 
ated by the desire to add adhesion into a Green's function mole- 
cular dynamics (GFMD) code used to model the contact 
between rough surfaces. To model adhesion, one needs to iden- 
tify a functional form for the finite-range surface forces 
allowing one to reach the JKR limit in an efficient and well- 
controlled fashion. It quickly became clear that doing a clean 
job is not entirely trivial and that modeling single-asperity 
contacts ought to be better understood first and moreover is 
interesting in its own right. 

The remainder of this paper is organized as follows: I first intro- 
duce the model, sketch the numerical methods and discuss diffi- 
culties arising in simulations in the limit of large and small 
Tabor coefficients. Next, I present a brief dimensional analysis 
motivating the commonly used unit system and the Tabor coef- 
ficient. The result section opens with adhesive contacts. There, I 
reproduce some established results and investigate how sensi- 
tive results are on the details of the interaction model. That 
section also contains a comparison to and an asymptotic 
analysis of the MD model motivating some minor modifica- 
tions of the empirical COS equations [1]. Next, results on posi- 
tive adhesion are presented before the major results are summa- 
rized and conclusions are drawn. 

Results and Discussion 
Definition of the model 

In this section, the single-asperity contact problem is intro- 
duced. As shown in Figure 1 , we consider a stiff, ideally-flat 
wall positioned in the xy plane at z = 0 and a linearly-elastic tip 
of parabolic shape. Its undeformed surface is given by 

z(r) = r 2 /lR, (1) 



where R is the radius of curvature and r = ^jx + y the 
in-plane distance of the center of the tip from the origin of the 
coordinate system. The elastic displacement of the tip, u(x,y) is 
formally a function of both in-plane coordinates, although the 
equilibrium solutions only depend on r. The gap g{x,y) indi- 
cates the distance between the deformed tip and the substrate, 
i.e., 

g{x,y) = z{r(x,y))-u(x,y). (2) 
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It is furthermore assumed that the tip cannot penetrate the sub- 
strate. This can be stated as a non-holonomic boundary condi- 
tion 

g{x,y)*o. (3) 

Alternatively, one can formally introduce a short-range, hard- 
wall repulsion [14] 

V SY = lim jd 2 r(f Y z v )exp{-g(x,y)/z Y }, (4) 

where f Y is an arbitrary positive constant of unit force per area. 
Note that the integrand on the r.h.s. of Equation 4 is zero for 
finite gaps while it diverges for negative gaps. Depending on 
the problem, it can be more convenient to use either the non- 
holonomic boundary condition or the energy-based description 
of the short-range repulsion. 



This work also considers finite-range adhesive or finite-range 
repulsive interactions Vf r , which only depend on the local gap. 
The default expression for it is: 

Vfr =-Ay{d 2 rexp{-^(x,7)/z 0 }, (5) 

where Ay is the work of adhesion per surface area that is 
obtained when a flat tip touches the substrate in a continuum 
description. The choice of the functional dependence of Vf r is 
not motivated by the true functional form for the interactions 



between any two real solids, but for the moment being, it is a 
matter of convenience. Alternative interaction models for the 
integrand are introduced in a seperate section. 

An important property of all models for Vf r is that the inter- 
action is characterized by a prefactor representing the work of 
adhesion and a single length scale zq. The choice of the latter 
allows one to localize adhesive stress near the contact line 
through zq — > 0 or to extend the adhesive interactions to radii 
much exceeding the contact radius a c for zq — > oo. Of course, zq 
can take any value between zero and infinity so that intermedi- 
ate-range interactions can be modeled as well. Note that V& and 
F sr are qualitatively different: The prefactor of the short-range 
potential is formally zero, because f r is finite and z r very small. 
In other words, z r represents the size of an "infinitesimally- 
small" atom whose size is irrelevant on a continuum scale. In 
contrast, the prefactor of the finite-range potential is considered 
finite as well as the range of interaction zq. It represents a 
"collective" length scale, such as the decay length of density 
oscillations in the fluid [15] or the radius of gyration of a 
polymer. 

The displacement u(x,y) and other fields (gap and stress) will be 
expressed not only in real space, but also in Fourier space. This 
is done most conveniently by using in-plane periodic boundary 
conditions. The respective boundaries lie atx or j; = ±L/2, which 
should be chosen such that L (the linear dimension of the simu- 
lation cell) is much greater than the linear dimension of the 
contact zone. The latter includes the contact and the area of 
non-negligible adhesive (or finite-range repulsive forces) 
stresses. The following convention for the Fourier transform 
shall be employed 

w(q) = — Jd 2 rw(r)exp{/q-r} (6) 

A 

w ( r ) = Z^) ex PW' r }> (7) 

q 

where the wave vector components satisfy q a = 2nn/L, A = L 2 is 
the integration domain, and n an integer number. With these 
definitions, one can express the elastic energy of the deformed 
tip (in the small-slope approximation) as 

F ela=^I>Kq)| 2 > (8) 

q 

where E is the effective modulus, E = Ey/(\ - v 2 ), Ey being the 
Young's modulus and v the Poisson ratio. The convention of 
using the symbol E* for the effective modulus is abandoned for 
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Figure 1 : Geometry of the deformed tip (upper grey solid), the sub- 
strate (lower solid), and the reference tip (solid line). The dotted line 
represents a hypothetical tip that is allowed to penetrate the substrate 
the distance d into the substrate without deforming. The following 
vectors are introduced in the sketch: Normal load Fn acting on the 
center of mass of the tip, the elastic displacement field u, and the dis- 
placement d of the tip's center of mass. In addition, two scalar quan- 
tities, namely the contact radius a c and the gap (field) g are shown. 
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clarity, because primes will be used later to indicate scaled 
coordinates and scaled parameters. 

Since t/(o) can be interpreted as the center of mass mode, the 
effect of a load (or normal force) exerted on the tip leads to an 
energy 



r«t=-*N«(0)- 



(9) 



When solving the contact problem, one seeks to minimize the 
total energy 



'tot 



[»(')] = 



^sr+^r + ^ela+^ext 



(10) 



with respect to u(x,y), i.e., the solution uo(r) must satisfy 



where ©(•••) denotes the Heavyside step functions, which is one 
for positive arguments and zero otherwise. The interaction 
potentials and their first derivatives are shown in Figure 2. 




5K 



tot 



[«(')] 



8w(r) 



= 0 



(11) 



\u(r) = u 0 (r) 



for each r. Here, 5 indicates a functional derivative. In a discrete 
representation of the problem, r is an index so that the func- 
tional derivative in Equation 1 1 has to be replaced by a partial 
derivative. 

Alternative interaction models 

Throughout this paper, different functional forms for the finite- 
range interactions between surfaces are considered in addition 
to the "default" or "exponential" model introduced in 
Equation 5. Functions similar to the ones used in this work have 
already been employed for the simulation of mode I fracture or 
delamination. Depending on the authors, the functions are 
called the cohesive zone model [16], the crack evolution func- 
tion [17], or the traction-separation relation [18]. However, it is 
not clear how the results obtained for mode I fracture geome- 
tries relate to Hertzian contacts. This is the main reason why the 
results obtained within the cohesive zone model cannot be 
compared in a straightforward fashion to those of the current 
study. 

The additional models in the current work replace the integrand 
on the r.h.s. of Equation 5 with the following expressions: 



-Ay x 



exp 
X 



(V/2,o 2 ) 
[l + UAof 



Gauss modell 



van der Waals modell , (12) 



(l-g/z 0 )©(z 0 -g) Maugis-Dugdale model 



Figure 2: (a) Finite-range surface energies and (b) forces per unit area 
for the models investigated in this study. 



All expressions take the same value, -Ay, for the adhesive 
energy when the two surfaces touch, i.e., in each case the work 
of adhesion is Ay. In this sense, all four models produce the 
same continuum limit. However, in two models, namely the 
Gauss and the van der Waals (VDW) integrands, the derivative 
v' fr (g/z 0 ) goes to zero when the two surfaces touch, while it 
remains finite in the exponential model and the MD model. 

As stated before, none of the models are supposed to be highly 
realistic representations of surface forces, although each model 
may have its justification. In particular the exponential model 
follows from the argument that long-range density correlations 
in fluids are governed by a single length [15]. In a high-density 
fluid, the correlation length becomes complex [15], which then 
leads to layering transitions as discussed recently [19]. The 
VDW model might approximate the long-range van der Waals 
interactions in a way that Ay reflects the Hamaker constant. 
Depending on the confined system in question, other effective 
interactions might be possible. However, all models represent 
the feature that surfaces repel upon close approach (i.e., when 
atoms from opposing surfaces bump into each other, which is 
implemented through the hard-wall repulsion) and that attrac- 
tion - or additional repulsion - may occur at finite distance. 
Continuous short-range repulsive forces are not used here. 
Doing so would complicate the definition of contact and thus 
contact radius, which has remained controversial for systems 
without hard-wall or hard-disk interactions [20-23]. Lastly, the 
equations to be solved would become very stiff and thus the 
simulation inefficient if the hard-wall constraint was replaced 
by potentials with large curvature. 
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In the context of the squeeze-out of fluids, the MD and the 
exponential model might not be physically meaningful for small 
ratios of g/zo: when one flat solid is placed on top of another flat 
solid with an infinite simally small external load (in the absence 
of a fluid), the two solids would repel, although they "cannot 
know", away from a contact line, that a fluid wants to penetrate. 
Nonetheless, the exponential model has been used in early study 
of squeeze-out of fluids [12]. Forces between two (flat) surfaces 
in the Gauss and the VDW model are zero either for intimate 
mecanical contact or at infinite separation. 

Dimensional analysis 

In this section, I present a simple dimensional analysis of the 
contact problem. The result of the analysis is a meaningful set 
of units, which, in similar form, has already been established by 
Maugis [9]. However, in the present analysis, units are not 
motivated from the solutions but rather straight from the begin- 
ning, i.e., by the expressions defining the model. This is why 
Maugis' and the present units differ by dimensionless prefac- 
tors, which, however, always turn out close to unity. In the 
subsequent derivation, it is not necessary to know the precise 
functional dependence of the finite-range forces. 

Assume we know the solution u$(y) minimizing V tot for a given 
set of parameters defining our model, i.e., u(r) solves the 
contact problem for a specific set of values for E, R, F^, Ay, 
and zq. It is then possible to "recycle" the shape of the function 
uo(r) to solve a different problem defined by different parame- 
ters E\ R', F{j, Ay', and z' 0 . Specifically, if each individual 
summand of ^t 0 t[ w '( r ')] i s identical to the equivalent term in 
F tot [w(r)] (up to a multiplicative constant, which can be set to 
one), then Wo( r ) minimizes ^t' 0 t[ w '( r )] given that wo( r ) mini- 
mizes V tot [u(r)]. 

The transformation, u' 0 (V) = tu(r'} = tu 0 (sr), in which in-plane 
coordinates are scaled as r' = sr and normal coordinates are 
scaled as z' = tz leaves the shape of the solution unchanged. Of 
course, z(r) and thus g(r) must be transformed the same way as 
w(r). Therefore, the radius of curvature of the "new" tip is 
R' = s 2 Rlt. 

Let us investigate how one has to alter each individual term 
contributing to F t ' ot [w'] so that it matches its equivalent in 
PtotM- (i) The hard-wall repulsive energy V& is unproblematic. 
It disappears for the old and the new solution, because of the 
limit z r — > 0, i.e., V[ x [w'] = 0. (ii) To recycle the V& calculation, 
we need to set z' 0 = tz$. Keeping in mind that A' = s 2 A, where 
A = L 2 is the original integration domain, it follows that 
= s 2 (Ay'/Ay)V[u]. (iii) For the calculation of the elastic 
energy, it is useful to keep in mind that q' = qls and that 
A' = s 2 A. This means that V^[u'] = st 2 (E'/E)V ela [u]. (The 



integer indices enumerating the wave vectors are identical for 
the original and the new domain.) (iv) Lastly, the load-related 
energy becomes F e ' xt =t{F^/F N "jV ext . In summary, we can 
recycle our solution with the following substitutions 



E' = E/st 2 (13) 

R' = s 2 R/t (14) 

Ay' = Ay/s 2 (15) 

z' 0 =tz 0 (16) 



Let us first consider the case of no external force, F^ = 0, so 
that we investigate the "intrinsic" system properties. If we use E 
as the unit of pressure, which is done until further notice, all 
three remaining parameters defining the system can be 
expressed to be of unit length, i.e, R, zq, and Ay IE. Wether a 
potential should be classified as short- or long-ranged can only 
depend on a non-dimensionalized interaction length. This 
means that zq has to be expressed in the two remaining units of 
length (R and Ay/E) such that the dimensionless ratio 

R l ' a (H/E) a (18) 

remains unchanged under a scaling transformation. 

Let us now chose the radius of curvature as the unit of length, 
or, alternatively, consider only those scaling transformations 
that leave R constant. This can be achieved by setting t = s 2 , 
which maps an infinite parabola (x — > z = x 2 ) onto itself 
(sx — > s 2 z). I note in passing that such a transformation might 
not be meaningful for a scaling analysis of the contact 
mechanics of randomly rough surfaces, which will be presented 
elsewhere. 

Inserting the relevant equalities from Equation 13-Equation 17 
reveals that choosing a = 2/3 leaves the expression in 
Equation 1 8 constant. As a consequence, the range of influence 
of the adhesive term is best quantified by the term 

^ = ~^^' (19) 
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where |ij is known as Tabor coefficient - up to a dimensionless, 
multiplicative constant. It remains invariant under all scaling 
transformation in Equation 13-Equation 17 leaving the radius 
of curvature unchanged. 

From Equation 19 one can see that we need to send zq — > w 2/3 zq 
in order to keep the Tabor coefficient constant when changing 
Ay IE at constant R to w Ay IE. This in turn implies a transforma- 
tion of x — > w 1/3 x for the in-plane coordinates, because R is 
supposed to remain unchanged. It follows that a c (F^ = 0) — > 
w 1/3 a c (F N = 0) and thus a c oc (lAyl/i?) 173 . The unit of a c can be 
fixed by multiplying the r.h.s. of this proportionality with R 2/3 . 
Otherwise the proportionality coefficient can only depend on 
and of course, on the sign of Ay. Therefore, we can write 



a c (F^=0) 



( r2 N/ e ) 



V3 



: Vt(^t^n=°)> 



(20) 



such that the r.h.s. of the equation only depends on the Tabor 
coefficient and the functional form of the surface interaction. 
Since we have not used the explicit functional form of our 
default surface interaction (other than that it depends only on a 
single length scale), the conclusions drawn in this section 
extends to any choice for V& considered in this work. 

To include finite loads into the analysis, note that the ratio 
F N /|Ay|i? does not change under the transformation 
Equation 13-Equation 17. This allows us to express a properly 
undimensionalized contact radius as a dimensionless function of 
a properly dimensionless load 



r/(R 2 Ay/E) 1/3 
2.1 2.2 2.3 2.4 2.5 2.6 2.7 




).08 0.085 0.09 0.095 0.1 0.105 0.11 

rls 



Figure 3: Scaled gap g(r) as a function of scaled distance from origin r 
for two different parameter realizations related through the scaling 
transformation Equation 13-Equation 17. Parameters used are 
E = L = u T = 1, F N = 1-10" 4 /wand Ay = 0.64-10~ 4 /w. The surface is 
discretized into 512 x 512 elements. Circles s = 1 with w= 1 — ► 
z 0 = 0.0016, crosses w = 1/2 z 0 = 2 2/3 ■ 0.0016. For the second data 
set, this implies a scaling factor of t = 2 2/3 for variables linear in normal 
coordinates and thus s = 2 1/3 for variables proportional to in-plane 
coordinates. The units on the normal side of the axes are in "absolute" 
units, i.e., L = 1 . The units on the opposite axes correspond to those 
that remain invariant under a scaling transformation. The dotted line is 
drawn to guide the eye. 



The normal displacement can be nondimensionalized in a 
similar fashion as the contact radius, except that it needs to be 
rescaled with the factor w 2/3 rather than with w 1/3 . This is why it 



must obey 



(22) 



ci c (F^ J / /i i \ 

1/3 =M^N/N4 (21) 

( R HI E ) 

From this last equation, it also becomes clear that the pull-off 
(or the squeeze-out) force is proportional to \Ay\R, i.e., by iden- 
tifying the value of F N /| Ay\R at which the function a c j takes its 
minimum value. Therefore, it is most meaningful to normalize 
the force with AyR, unless, of course, Ay = 0. 

The approach is validated in Figure 3. It shows the spatial 
dependence of the gap for two different parameterizations. 
Small deviations, which are not visible to the naked eye, occur. 
They are due to finite-size and discretization effects. For 
example, the ratio aJL is not exactly zero and takes different 
values for different values of s for a fixed number of points used 
to represent the elastic surface. 



where all terms on the r.h.s. of the equation are again dimen- 
sionless, while those on the l.h.s. are allowed to have units. 
Thus, displacements and gaps are best represented in units of 
i? 1/3 (|Ay|/£ T ) 2/3 , while contact radii are more meaningfully 
expressed as multiples of i? 2/3 (|Ay|/is) 1/3 . As a result, numbers 
turn out of order unity when data is represented in these units, 
unless F N approaches the pull-off threshold or distinctly 
exceeds AyR. 

I conclude this section by summarizing the units used in this 
work and discuss some of the consequences arising from it: 
Specifically, the following units are used for: 

in-plane coordinates: [x] = R 2 ^ {^Ay\l (23) 

out-of-plane coordinates: [z] = R l/3 (|Ay|/^) 2/3 (24) 
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normal force: [^n] = ^|^y| (25) 

normal pressure: [p]= (|Ay|/tf) 1/3 £ 2/3 (26) 

This list includes a "new" unit of normal stress or pressure, [p], 
which must be chosen proportional to [FnMx] 2 rather than to E 
so that the regular definition of pressure applies. As noted 
above, E drops out of the definition of the unit for the normal 
force, implying that pull-off or squeeze-out forces cannot be 
functions of E. Instead they must equal R\Ay\ times dimension- 
less expressions that can only depend on jij, and, of course, on 
the functional form of the interaction potential. In our units, the 
well-known a c (F^,E* ,R,Ay) relations can be simplified to 

a\ = — (F N + In) in DMT limit (27) 



a\ =-[V N +37i + ^67tF N +97I 2 J in JKR limit. (28) 

Hertzian contact mechanics is obtained in either limit for 
» 1. Finally, note that Maugis' choice for units slightly 
differs from ours in that he used rcAy rather than Ay in 
Equation 23-Equation 26 and 3E/4 instead of E. The conver- 
sion between Maugis' and our system is summarized in 
Equation 43-Equation 46. 

Numerical analysis 

Different methods can be used in the numerical solution of 
Equation 11. For the present study, Green's function molecular 
dynamics (GFMD) [24] as described in [25] was employed. The 
only modification is the implementation of conservative surface 
forces acting in addition to the boundary condition g(r) > 0. 
Moreover, the results in this work were produced with a serial 
code with typical run times of a few minutes. I refer to the 
literature for more details on GFMD [24,25]. Irrespective of the 
employed code or method, particular precautions, which are 
worth discussing, need to be taken into account when including 
adhesion or finite-range repulsion. 

When simulating Hertzian contact mechanics, one needs to 
ensure that the discretization of the lattice Aa satisfies Aa « a c . 
Of course, methods based on spatially varying grids only need 
to obey that relation near the contact line. In addition, one wants 
a c to be much less than the size of the simulation cell, at least in 
Fourier-based methods, such as GFMD. One then has the 
sequence of inequalities Aa « a c « L. In Hertzian contact 
mechanics, this is easy to achieve: choosing the discretization 
such that Aa/a c = 1/32 and a c IL = 1/8 already gives accurate 



results for the contact area, that is, to within less than 0.1% 
error if the contact area is determined through a fit of the radial 
pressure profile. 

When including adhesion, an additional length enters, namely 
that associated with the adhesive zone. The additional adhesive 
radius or skin a a then needs to be taken into consideration. 
When the Tabor coefficient is very small, a a becomes large, and 
one needs a a to lie within the simulation cell. A new series of 
inequalities is obtained: Aa « a c « a a « L. If the normal 
stress changes smoothly with the gap, i.e., for long-range adhe- 
sion, the forces couple predominantly to large wavelength 
modes. This then results in a simple offset force, as is well 
known from DMT theory. However, numerical demands can 
become significant when a c disappears continuously with 
decreasing load as in the DMT scenario. The condition a c « a a 
then becomes difficult to satisfy. 

In the opposite case of a large Tabor coefficient, a a is very 
small, potentially much smaller than a c . We still need to resolve 
the adhesive zone sufficiently well, because the stress has to be 
smooth on the given discretization. One thus obtains the series 
of inequalities a « a a « a c « L. In either limit of large or 
small Tabor coefficient, another inequality needs to be satisfied 
in addition to those for Hertzian contact mechanics. 

While the contact area converges reasonably quickly as the 
respective inequalities are obeyed, the center-of-mass displace- 
ment d, which corresponds to w(o) or u(r — > oo) only converges 
slowly. The reason is that the displacement field induced at a 
given point due to an external force only decays with the 
inverse distance from that point, i.e., 

d&u(r) + — for r»a c +a a , (29) 
r 

where c is a load- and system-dependent constant. Outside of 
the adhesive zone, this relation can be used, in principle, to 
extrapolate from finite r to infinite r, i.e., by determining c and 
d from two measurements taken sufficiently far outside the 
adhesive zone. In practice, this turns out problematic, because 
the periodic boundary conditions suppress the \lr corrections 
near the boundaries. 

For Hertzian contacts, it is still possible to use a slightly modi- 
fied (empirical) correction 

d&6u x -5u M , (30) 

where X and M denote the (symmetry) points (L/2,L/2) and 
(L/2,0) relative to the position of the center of the tip. The same 
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extrapolation scheme also appears to give quickly converging 
estimates for the normal displacement for adhesive contact, 
which is demonstrated in Figure 4. 



It is worth discussing Figure 4. At the given load, the contact 
radius is a c ~ 2.30, while it would have been a c ~ 1.05 without 
adhesion. The displacement curve has a peak at r = 2.52 and 
adhesive effects remain non-negligible all the way up to r ~ 4. 
At that distance the gap starts to be greater than 5zq, which 
means that the adhesive attraction is less than e~ 5 times the 
value in the contact and its immediate periphery. For distances 
exceeding r = 4, an infinitely large system would then show the 
displacement given in Equation 29. The periodic boundary 
conditions suppress this scaling rather strongly, yet, for radii as 
small as r = 10, accurate estimates for can be achieved 
through Equation 30. 

Simulations could be made more efficiently by exploiting the 
radial symmetry of the system. This would allow one to reduce 
sums over two indices (e.g., q\ and qj) to that over one index. 
However, less is gained than it first might seem. To get a good 
resolution of the contact area, the one-dimensional (ID) calcu- 
lation require greater ratios for aja than two-dimensional 
setups. The reason is that the resolution of the contact area in 
ID and in 2D both scale with l/N, where N is the number of 
grid points in the contact. For example, when representing a 
contact in which for the given discretization 5a < a c < 6a in a 
2D system, then a c is allowed to take the values V5 2 +l 2 , 
V5 2 +2 2 , V5 2 +3 2 , and V4 2 +4 2 . The maximum distance 



between two radii thus is Aa max = 0.28 so that the resolution is 
Aa c /a ~ 5/0.28. To match this in a ID model, one would need 
18 grid points rather than five. Since we need to cover a 
(square) area of (2a c ) 2 in 2D, we thus have a computational 
overhead of a little more than a factor of 4 compared to ID. 
However, the disadvantage of large ID systems is that the 
number of iterations to find solutions can be much larger than in 
2D. Depending on the nature of the solver, the number of itera- 
tions scales as a power law with linear system size. In the given 
case, where the effective stiffness scales proportional to wave 
vector q, one would expect a slowing down with yfq, at least in 
simple gradient-based minimization such as steepest descent or 
molecular-dynamics. For this reason, no efforts were made to 
reduce the dimensionality, although this would have been legiti- 
mate for the given problem. 

Positive work of adhesion 

This section analyzes how the employed models reproduce 
established results for adhesive single-asperity contacts in the 
limits of large and small Tabor coefficient. This includes an 
asymptotic analysis of Maugis' solution of the Dugdale model, 
which in turn leads to modifications of the equations proposed 
by Carpick, Ogletree, and Salmeron [1]. The cross-over from 
JKR to DMT is investigated as well, in particular at zero load 
and near pull-off, allowing one to identify the model for the 
surface interaction that is most appropriate for the simulation of 
(adhesive) multi-asperity contacts. 

Zero external load 

An external load of = 0 is addressed first. The motivation for 
studying this special case is that one can analyze relatively 
easily at what Tabor coefficients the DMT and JKR limits start 
to predict reasonably accurate values for the contact radii and 
displacements in our various models. 

We start our analysis with the pressure distribution of the expo- 
nential model, which is depicted in Figure 5 for jij = 1/4 and 
jix = 4. It behaves very similar to the MD model, which shall 
not be shown explicitly. As to be expected, the adhesive load is 
spread out for |i T = 1/4 to radii clearly exceeding a c (all the 
more as each radius r contributes with a weight proportional to 
r), while it is rather localized near r = a c for = 4. It therefore 
seems legitimate to call the (net) pressure profile for jij = 1/4 
DMT-like and JKR-like for |i T = 4. 

The adhesive pressure is calculated from the functional deriva- 
tive p a dh( x >y) = ~SFf r /5g(x,3;), where Vf Y is defined in 
Equation 5. This can be evaluated to yield 

^adh^) = — ex P{-#( r )/ z o}> (31) 
z 0 




4 6 8 10 12 14 16 



Figure 4: Negative elastic displacement -u(r) (as defined in Figure 1) 
of a tip with R = E = 1 pressed against a rigid substrate at a positive 
load of Fn = 1 -5625 for a work of adhesion Ay = 1 and a Tabor coeffi- 
cient uj = 1 resulting in an exponential decay length of zq = 1. In each 
case, the system is discretized into 512 x 512 grid points, but different 
sizes are used, i.e., L = L 0 ,L = 2 V3 L 0 , L = 2 2/3 L 0 with L 0 = 9.8825. Part 
(a) shows a larger domain including the shape of the tip in form of a 
dashed line. Part (b) shows a smaller domain and includes a higher- 
resolution estimate of the displacement at infinite radius R through the 
extrapolation 6l/x - 5i/m- 
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Figure 5: Interfacial pressure on a free, linearly-elastic half space 
resting at zero external load on a rigid, adhering, parabolic substrate 
for (a) uj = 1/4 and (b) uj = 4 in case of the exponential model. In 
each case, the upper and lower grey lines indicate, respectively, the 
pressure due to the constraint and that due to adhesion. Contact radii 
are indicated by a c . The net pressure is shown by a black line as well 
as by small circles representing the actual grid points. Units are 
chosen according to Equation 23-Equation 26. 



which becomes p a dh( r < a c) = _ Ay/zo within the true contact 
area. Using Equation 19, one obtains 



R 



in physical units 



(32) 



-Mr 



in our reduced units. 



Stress or pressure originating from the constraint g(x,y) > 0 is 
computed from the elastic Green's functions. 

The well-known qualitative difference for the contact geometry 
of systems with large (jix = 4) and small (|ij = 1/4) Tabor coef- 
ficient is borne out in the radial dependence of the gap g(r). 
Specifically, Figure 6 reveals that a small Tabor coefficient 
makes g(r) have a positive curvature at r > a c , as in the DMT 
solution, while it has a negative curvature at r > a c , indicative of 
an adhesive neck, for large jij- Figure 6 also shows that the dis- 
placement (defined by the difference between the actual gap and 
the gap in an undeformed contact at r » 1) is smaller for 
jix = 4 than for jij = 1/4, although the contact radius is larger 
for u. T = 4 than for |i T = 1/4. 

The Gauss model behaves qualitatively different from the expo- 
nential model. First, there are no adhesive forces within the 
contact, but only outside of it, as shown in Figure 7. Second, at 
r = a c , the total pressure disappears in the Gauss model, while it 
remains finite in the exponential model. Third, the pressure due 
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Figure 6: Gap between a rigid, adhesive, parabolic tip and a linearly- 
elastic half space for two different Tabor coefficients uj = 1/4 (solid 
black line) and uj = 4 (broken black line). The gap z = x 2 I2 of an unde- 
formed half-space is shown in grey for comparison. Arrows indicate 
contact radii. Units are chosen according to Equation 23-Equation 26. 



to the constraint has finite slope at r = a~ in the Gauss model, 
while the slope diverges in the exponential model (not shown 
explicitly). All these differences arise because the derivative of 
Vf r (g/^o) remains finite (i.e., with positive sign) in the exponen- 
tial model when the gap closes, while Vg.(g/z 0 ) is zero in the 
Gauss model. 

Another consequence of p a dh(g/ z o) having zero slope in the 
limit of g — > 0 + is that the gap in the Gauss model closes con- 
tinuously rather than with a discontinuity in its first derivative. 
This is shown in Figure 8, from where it becomes clear that it is 




Figure 7: Pressure p(r) in the Gauss model at zero load for two 
different Tabor coefficients as a function of the in-plane distance rfrom 
the center of the contact. The contributions due to the constraint are 
positive, i.e., above the grey line, while the adhesive forces are nega- 
tive. 
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difficult to define good measures for the contact radius. In a 
linear representation (at low resolution), it seems as if the 
contact closes with the typical adhesive neck, i.e., in part (a) of 
Figure 8 the gap appears to close at r ~ 2.415. There, the slope 
of g(r) takes its maximum value, in a very similar fashion as in 
the JKR limit, or for the exponential model for the same value 
of [ij= 16. However, when increasing the magnification, one 
can see that the contact closes only at r ~ 1.97. Unfortunately, 
the radii where the gap closes to zero, and where g(r) has its 
maximum do not approach each other quickly when |i T is 
increased. Instead, the value of g in the cross-over region in 
Figure 8b moves to smaller values as |ij increases. Similar 
behavior is seen in the VDW model, which is not shown here 
explicitly. 



10 

K 

10 
10 

16 
12 



. 1 1 1 1 1 1 1 1 1 

: (a) 1 


1,1,1,1,1,- 


- , i , i , i , 


i 




1 




I 1 i 1 i 1 i 1 
I (b) 

1- i i i i , i i 


ri 


1 1 1 ! 1 ' h — 1 1 ! 

— Gaussian model -j 
■-■ exponential model -j 

I , I i I , I i I i -1 



. 1 1 1 1 1 1 1 1 1 
[ (c) 

i i i i i i i i i 


i 1 i 1 i 1 i 1 i 1 . 

L___ -i 


i , i i i , i i i i 



0 0.5 1 



1.5 



2.5 
r 



3.5 



4.5 



Figure 8: Gap g(r) in linear (a) and logarithmic (b) representation as 
well as (c) first derivative g\r) for p T = 16 in the Gauss model (straight 
lines). A higher resolution representation of the gap is given in (a) with 
a dotted line. The exponential model is shown for comparison in (b). 
The two thin grey lines are drawn where the gap becomes zero (left 
line) and where the slope of g(r) takes its maximum value. 



Zero-load contact radii for different potentials are depicted in 
Figure 9 as function of the Tabor coefficient. In the exponential 
model, the contact radius approaches DMT and JKR limits in a 
very similar fashion as in the MD model. In a later section on 
the asymptotic analysis, I find that the MD corrections to the 
JKR limit are of order l/ for large Tabor coefficients while 
those to the DMT limit are of order |i T for small Tabor coeffi- 
cients. The same scaling of the leading-order corrections is 
apparently borne out in the exponential model. 

Models in which Vf r (g/zo) has zero slope in the origin behave 
qualitatively different from the MD or the exponential model. 
They approach the DMT limit for the contact radius fairly 



— Maugis solution 
• exponential model 
x Gauss model 
□ VDW model 




2.5- 
2.4- 
2.3- 
2.2- 

2 2A \ 

s° 2.0 - 
1.9- 
1.8- 
1.7 
1.6 



Figure 9: Contact radius a c at zero external load as a function of the 
Tabor coefficient uj for different model interactions: exponential (full 
circles), Maugis Dugdale (straight lines), Gauss (crosses), and VDW 
(open squares). The upper and lower broken line denote the DMT and 
the JKR limits, respectively. In the Gauss and VDW models, two 
different definitions are used for a c : The upper symbols refer to the 
position where g\r) reaches a local maximum, while the lower symbols 
indicate the points of first gap closure (g = 0). 



quickly, i.e., roughly with ji^ 2 . However, convergence to the 
JKR limit is poor. The latter can be improved by defining the 
contact line to be located where g'(r) takes its maximum value. 
Unfortunately, this definition cannot be universally applied, i.e., 
only when |i T is sufficiently large to allow for an adhesive neck 
to be formed, see also Figure 8. Moreover, in the context of 
randomly rough surfaces with complicated contact geometries, 
this last definition of contact would not be practicable. 

Unlike the contact radius, the normal displacement d does not 
suffer from any difficulties to be properly defined. In principle, 
this could enable one to ascertain Vf r (g) from displacement 
measurements without much ambiguity. However, Figure 10 
reveals that the functional form of <i(ur) is relatively insensi- 
tive to the details of the finite-range interaction, at least, as long 
as we allow for a redefinition of the Tabor coefficient, such that 
all curves superimpose at the distance half way between the 
JKR and the DMT limit. This is in agreement with a work by 
Tvergaard and Hutchinson [18] who found that Ay and the peak 
stress (which one may losely associate with Aj/zq) are the basic 
parameters for mode I fracture. 

Before proceeding to the case of finite load, I wish to comment 
on the relatively large numerical (GFMD) uncertainties for the 
displacement at large Tabor coefficients. They stem predomi- 
nantly from the difficulty to apply the finite- size extrapolation 
formula, Equation 30, to gaps having adhesive necks. This 
problem would not be present in large-scale simulations of 
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Figure 10: Normal displacement of at zero external load as a function 
of the Tabor coefficient uj for different model interactions: exponential 
(full circles), Maugis Dugdale (straight lines), Gauss (crosses), and 
VDW (open squares). The Tabor coefficient was normalized so that all 
curves superimpose at the "half-way" distance d = (c/jkr + dowj)/2. 
The upper and lower broken line denote the DMT and the JKR limits, 
respectively. Numerical uncertainties arise in the limit of large Tabor 
coefficients, as indicated by the error bar. 

multi-asperity interfaces, because system sizes would automati- 
cally be much larger than local contact radii. One may conclude 
that the use of the exponential model for the study of adhesive 
multi-asperity contacts appears to be appropriate. The MD 
model could be used as well, in principle, however, it might in- 
duce undesired numerical instabilities due to the discontinuity 
of v'iglzo) at g = zq. The Gauss model can only be taken when 
the property of interest is related to the gap but not for the 
calculation of contact area. If one wanted to simulate van der 
Waals attraction at large one might want to replace the 
VDW model in Equation 12 with 1/(1 + g/zn) 2 . This depend- 
ence makes it possible to determine the contact area meaning- 
fully in the realm of continuum mechanics while using reason- 
able approximations for van der Waals interactions at large dis- 
tance. 

Finite external load 

In most experiments, the Tabor coefficient is kept constant and 
the normal load is changed. As a result, one obtains the normal 
displacement d(F^) as a function of the normal load F^. In 
some cases, i.e., for sufficiently large contact radii, an estimate 
of the contact radius, a c (F^), can be obtained as well. One 
might be tempted to believe that knowing such curves allows 
one to deduce the surface forces. Here, I want to investigate to 
what degree such an inversion is possible by studying the sensi- 
tivity of the functions d(F^) and a c (F^) to the functional form 
of the surface interactions. Figure 1 1 shows the contact radius 
as a function of the normal load. One can see that the exponen- 



tial model and the MD model agree very closely, that is, curves 
almost superimpose for a given Tabor coefficient. This makes it 
essentially impossible to discriminate between these two forms 
of interaction experimentally. Likewise, the Gauss and VDW 
models also coincide for the same Tabor coefficient despite 
their significant differences at large gaps. Interestingly, the 
|ij = 1 curve for VDW and Gauss (both having finite slope 
potentials at g = 0) is akin of the [ij = 1/4 curves for the MD 
and the exponential model (both having zero- slope potentials at 
g = 0). This confirms the trend reflected in Figure 9: Surface 
potentials with zero slope at g = 0 make the results move toward 
the DMT limit. 

In contrast to the a c (F^) dependence, the normal displacement 
curve <i(F N ) predominantly depends on the Tabor coefficient. 
Now, all jix = 1 curves resemble each other closely, inde- 
pendent of the slope of the surface potential at zero gap. As for 
the normal displacement, all curves are reasonably close to the 
JKR limit. Even the |ij = 1/4 curve lies closer to the JKR than 
to the DMT line. This is consistent with the results shown in 
Figure 10, which show that the DMT limit for d(F^) is only 
reached at extremely small Tabor coefficients. Figure 12 reveals 
that it is possible to adjust the free parameters of the MD model 
to fit d(F^) curves for a broad variety of surface interactions. 
However, one should abstain from deducing contact areas based 
on such fits, as this can result in non-negligible errors. For 
example, if we only knew the contact area from Maugis' solu- 
tion, we would be ill-advised to conclude from Figure 12 that 
the contact area for the jij = 1 Gauss model should lie roughly 
half way between those of the |ij = 1 and [ij = 4. 
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Figure 1 1 : Contact radius a c as a function of load Fn f° r the exponen- 
tial and the MD model using different Tabor coefficients, ranging from 
u T -> 00 (JKR, top) to u T = 0 (DMT, bottom). For the Gauss and VDW 
models, only uj = 1 is shown. Their a c (F N ) curve is similar to that of 
the Maugis and the exponential model for uj = 1/4. Color coding: 
uj = 4 (red), uj = 1 (green), and uj = 1/4 (blue). 
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g 



(m) = \jm 2 -\ acos(l/m)-m + 1 (38) 
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Figure 12: Normal displacement d as a function of load F|_ for the 
exponential and the MD model using different Tabor coefficients, 
ranging from uj = 0 (DMT, top) to uj —*■ 00 (JKR, bottom). For the 
Gauss and the exponential model, data is only shown for uj = 1 ■ Color 
coding: uj = 4 (red), uj = 1 (green), and uj = 1/4 (blue). 



and 



m 2 -l. 



(42) 



Comparison to other models and asymptotic 
analysis 



In each but one (straightforward) case, behavior of the func- 
tions for m approaching unity or infinity has been included. 
Maugis proposed an analytical solution for the relation between They become useful in the limit of large and small Tabor coeffi- 
contact radius a c and normal load in the Dugdale model. It cients, respectively, 
requires the elimination of an auxiliary variable, m, through the 

self-consistent solution of two coupled non-linear equations. Conversion back to our unit system can be done using: 
Once a c and m are found, the displacement d can be readily 
calculated as well. Using tildes to indicate variables in Maugis' 
unit system, the relevant equations read: 



1 = Y/W + 3 



(33) 



(34) 



F N = 5 3 - ja T a 2 /z(m) 

d = a 2 - d[i T j : {m) , (35) 
where the functions j\m), g(m), h(m), and j(m) are defined as 
f(m) = yjm 2 -I + (m 2 - 2^acos(l/m) (36) 



Mt = {2tt/9} 1/3 ft T 
a = (3tt/4) 1/3 a 
d = (3n/4) 2/3 d 



(43) 

(44) 

(45) 
(46) 



To overcome the need of having to find the self-consistent solu- 
tion to Maugis' equations, Carpick, Ogletree, and Salmeron 
(COS) [1] proposed a simple and thus elegant analytical 
formula for the a c (F^) dependence 



2(m-\)y formal 
— ^m 2 - 1 j form ^oo 



a c (0,n T ) 



2/3 



l + a(n T ) 



. (47) 



(37) 



Schwarz [2] later recognized that the COS description is exact - 
given proper parameterization - if the interaction between the 
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surfaces results from the superposition of an infinitesimally 
short-ranged and an infinitely long-ranged contribution. 
However, in the given context of one intermediate-range poten- 
tial, I will treat the COS equation as a guessed approximation 
containing the correct functional form in the limits of large and 
small |ix- 

The primary COS equation (Equation 47) is designed such that 
the contact radius at zero load a c (0,[ij) as well as the pull-off 
force Fpdij) can be reproduced exactly. However, approxima- 
tions to their dependence on ji T had been provided as well, 
because no closed-form expression are available. A free para- 
meter remains, a(|ix), which can be used to minimize devia- 
tions from the exact solution. At large loads, one recovers the 
well-known a c oc scaling, however, not necessarily with 
the correct prefactor. Another property of the COS approxima- 
tion is that it does not necessarily contain the correct value of 
the contact radius at pull-off. Thus, despite predicting the 
contact radius fairly well, the COS contact radius is not exact in 
the limit of very large and very small (i.e., pull-off) normal 
loads. These deficiencies can be improved when parameter- 
izing the COS equation in a slightly different fashion, e.g., 

a c 3 (F N ^ T ) = |(F N +F 0 (^ T ) + a( MT ) > /F N +F p ( MT )) (48) 
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Figure 13: Relative errors in per cent for the contact radius a c for 
ut = 1 at (a) negative and (b) positive load. The full line indicates the 
error when using Equation 48, while the dot-dashed line is based on 
Equation 47. In both cases, the parameter a(uj ) was adjusted to mini- 
mize the deviation from Maugis solution in the domain -F p < Fn < 2F p . 



tive errors are small, i.e., <1% even for [ij ~ 1, where one is 
relatively far away from both the DMT and the JKR limit. 

Zero load: The asymptotic analysis is readily done for zero 
loads, because the variable m can be directly eliminated in that 
case. As a result, one obtains 



with 



4 ^ 

F 0 (^ T ) = F p (|i x ) + -a c (0,^ x ). 



(49) 



5 c(^N=°) = 



6 1 / 3 + 0^i T 2 j for |li t co 
2 1/3 +0(|u T ) for|I T ^0 



(50) 



This set of equation ensures that a c converges to the exact value 
when F N — > oo and F N — > -F p . The parameter a(|ix) can then be 
adjusted to either yield the correct zero-load contact radius, or 
to minimize the deviation between approximation and the exact 
Maugis solution by some other mean. Note that the factors 3/4 
in Equation 48 and 4/3 in Equation 49 have to be replaced by 
unity when working with Maugis' unit system. 

I wish to note that including the correct asymptotics in the 
a c (F^) expression does not necessarily improve the fits in the 
range from slightly above the pull-off force at negative loads to 
several times the absolute pull-off force. This is demonstrated in 
Figure 13. Moreover, convergence to the correct a c (F N ) 
dependence at large loads is rather slow even when using 
Equation 49. 

It is also possible to constrain the COS relation for the contact 
radius such that it contains the correct pull-off force and contact 
radius as well the correct zero-load radius. In either case, rela- 



and 



— » OO 



(4/3) 1/3 +C»(fii 2 )forji T 
2 2/3 +C>(fi^ 2 ) forji T ^0 



(51) 



From the last two equations, one can see - as in Figure 9 and 
Figure 1 0 - that the JKR limit is quickly reached as the Tabor 
coefficient increases. However, convergence to the DMT limit 
with decreasing ju T is rather slow. It is particularly slow for 
the normal displacement. E.g., to have a maximum error in 
a c (F N = 0) and <i(F N = 0) that is of order 1 % with respect to a 
desired limit, it is sufficient to work with jLx x ~ 10 for JKR, but 
one needs ju x ~ 1 0~ 4 to approach the DMT limit with similar 
accuracy. The latter is not problematic for the simulation of 
multi-asperity contacts, as the system is large by default. 
However, for single-asperity contacts, large deviations from 
|i T = 1 (on a logarithmic scale) are difficult to handle in single- 
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asperity contact simulations for reasons discussed in the nu- 
merical-analysis section. 

Knowing the asymptotic behavior of 5 0 (jli t ) = 5^F N =0, jI T ) 
and d 0 (fi T ) = d (f n = 0, jii T ) with respect to ju x allows one to 
incorporate it in empirical equations for these two quantities. 
The following equations are found to achieve this and to 
provide excellent approximations to Maugis' solutions: 



1/MT +c 4 +c 3Mt 



(52) 



( d DMj/ V^T + c 5 )/^T 3 + C 7 ( c 6 + ^JKR^T ) 
d 0 (jI T J~ — (53) 



|Li T 5/6 + Cl \i\ 



Two coefficients in each of the last two equations (c\, c 2 and C5, 
eg) can be constrained to reproduce the correct asymptotics (and 
thus be obtained analytically). Two fit parameters remain for 
contact radius and one for the displacement. The relative errors 
from the pertinent fits are shown in Figure 14. Compared to an 
already quite accurate empirical relation proposed by Carpick et 
al. for a 0 (jLi T ), see Eq. (12b) in [1], the new Equation 52 and 
Equation 53 contain the correct asymptotics and reduces the 
maximum relative error from 1.5% to 0.3%. The data shown in 
Figure 14 were obtained with the following numerical values: 
c x = 4/5, c 2 = -1.285, c 3 = 4/5, c 4 = -0.435, c 5 = -3/2, 
c 6 = 0.1845, and c 7 = 6.71. 
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Figure 14: Relative errors in per cent for (a) contact radius and (b) 
normal displacement at zero normal load. Full lines refer to a fit based 
on Equation 52 and Equation 53 containing the correct asymptotic 
limits. The dotted line reflects an empirical fit based on the COS equa- 
tions. 



Asymptotic behavior near pull off: The structure of the COS 
approximation, Equation 47, and its modified form, 
Equation 48, indicates that the critical behavior near pull off 
satisfies a c - a p oc + F p j , where k must be either k = 1/3 
as in the DMT limit, or k = 1/2 as in the JKR limit. However, 
nothing in the self-consistent solution of Maugis indicates that 
the exponent k changes discontinuously from one value to the 
next as the Tabor coefficient reaches or passes through a crit- 
ical value. In fact, representing the data from Figure 1 1 in terms 
of A« c (AFn), as done in Figure 15, shows that k changes con- 
tinuously from 1/3 to 1/2 as \ij increases from 0 to infinity. 
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Maugis solution )i T = 1/4, 1, 4 
exponential model 
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Figure 15: Excess contact radius Aa c = a c - a p as a function of the 
excess load AFn = + for different values of the Tabor coefficient 
ranging from uj = 0 (DMT, top) to uj — > 00 (JKR, bottom). Here a c and 
a p denote the contact radius at an arbitrary load Fn and the pull-off 
load F p , respectively. Deviations between the Maugis solution and the 
exponential model are particularly obvious for the uj = 1 data set. 
Color coding: uj = 4 (red), uj = 1 (green), and uj = 1/4 (blue). 



An analysis for the normal displacement, similar to the one 
presented in Figure 15 but not shown explicitly here, exhibits a 
similar trend. The exponent describing Ad = d - <i p as a func- 
tion of AFn = + crosses over continuously from the DMT 
to the JKR limit as \ij increases. However, there is not a one-to- 
one relation between jij and k. In particular the data sets for 
jix = 1 show relatively large differences between the exponent 
in the MD model (k ~ 0.469) and the exponential model 
(k ~ 0.435). 

The insights obtained from Figure 1 5 can be used, in principle, 
to further modify the COS approximations, e.g., by replacing 
the square-root in Equation 48 by some other power or likewise 
by changing the square-root and the exponent 2/3 on the r.h.s. 
of Equation 47 in an appropriate fashion. When doing so, the 
modified version of Equation 48 does not only converge to the 
correct value at pull-off. It can also be parameterized to yield 
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the correct asymptotics near pull-off. This results in a further 
reduction of the mean or overall error by a little more than a 
factor of two with respect to those shown in Figure 13, 
however, at the expense of one additional fit variable. Since the 
main new aspect of this study is concerned with negative work 
of adhesion and, moreover, both original and modified COS 
equations are already quite accurate, a more detailed analysis of 
the adhesive single-asperity contact is not pursued in this work. 

Negative work of adhesion 

For repulsive contacts, Ay < 0, there is obviously no finite 
contact radius at zero normal load = 0 + . The repelled rigid 
tip simply "hovers" at (infinitely) large distance over an unde- 
formed elastic manifold. This is why it is not possible in this 
case to conduct a zero-load analysis similar to that presented for 
adhesive contacts. Since Maugis' solution has not yet been 
extended to repulsive contacts, we are not in a position to 
compare our data to analytical solutions for negative Ay. One of 
the consequences is that the asymptotic analysis must be based 
on GFMD data, except for [ij — > 0, for which normal forces 
couple predominantly to long-wavelength modes so that the 
Hertz-plus-offset approximation (DMT) should be accurate. 
Given the close similarity between the exponential and the 
Maugis-Dugdale model as well as that between the Gauss and 
the VDW model seen in the last section, the attention is 
restricted to one potential in each class, i.e., the exponential and 
the Gauss model. 

We start our analysis with the contact radius dependence on 
load. In analogy to the context of wetting fluids, one may call 
the force at which a finite value of a c becomes unstable upon 
lowering the load the spontaneous wetting force F sw . The force 
above which a c can no longer be zero is called the squeeze-out 
force F sq . If the transition from contact to non-contact is con- 
tinuous F sw = F sq , otherwise F sw < F sq . Results are shown in 
Figure 16. 

As is the case for attractive interactions, the contact radius at 
small loads can be sensitive to both the Tabor coefficient and 
the choice of the potential. Specifically, the exponential model 
always shows a continuous transition from finite to zero contact 
radius (at least for the values of |ij investigated here), while the 
Gauss model has either a continuous transition below a critical 

Tabor coefficient |u T < 1 or a discontinuous transition for 

* 

ji T > ju T . The discontinuity of the contact radius for Gauss 
potentials and sufficiently large Tabor coefficients implies that 
two solutions may coexist, i.e., one where the two surfaces are 
separated and one where they touch. However, once exceeds 
a second threshold force F S q(\ij), i.e., the squeeze-out force, 
only one solution survives, that is, the one with finite contact 
radius. This can be seen in analogy to adhesive contacts with 




Figure 16: Contact radius a c as a function of normal force F N for the 
exponential (full symbols) and the Gauss (open symbols) model. Lines 
connect data points (not all shown explicitly). In the case of uj = 4 
(Gauss model), an arrow indicates where the a c = 0 solution becomes 
unstable for increasing F N . Color coding: uj = 4 (red), uj = 1 (green), 
and uj = 1/4 (blue). 

jix > 0, where two solutions coexist in a finite interval of forces 
-F p <F N <0. 

As for the a c (F^) relation near pull-off in the case of positive 
work of adhesion, the excess contact radius, a c - a sw , depends 
as a power law on the excess force, F^ - F sw , for F^ > F sw : 

*c-*sw x ( F N- F Sw) K - ( 54 ) 

Fits to the a c (F^) relation are shown exemplarily for two values 
of |ij in Figure 17. Details about the fits to the presented as well 
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Figure 17: Contact radius a c as a function of normal force F N ' n the 
vicinity of the spontaneous wetting force F sw . Symbols reflect numeri- 
cal results while lines are fits according to Equation 54. They termi- 
nate at a c (F sw ). In the case of uj = 1 , an arrow indicates where the 
a c = 0 solution becomes unstable with increasing F N . 
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as additional data are summarized in Table 1 . As for attractive 
contacts, it is found that k changes continuously from 
K (HT ~~ > °°) = 1/2 to k(|1x —> 0) = 1/3. For small jij , Hertz-plus- 
offset behavior is reached as evidenced by the observation that c 
and F sw approach (3/4) 1/3 and 2n, respectively. However, F sw 
as well as F§q quickly increase with \ij for jix > 1 . This latter 
behavior is different from that of the pull-off F p force for attrac- 
tive surfaces, which only varies between 1.571 and 2n in the 
present unit system. Since the increase of both F sw and F sq is 
much faster in the exponential model than in the Gauss model, 
one can conclude that the exponential model converges more 
quickly to the continuum model than the Gauss model. 



Table 1: Results of fits to the data shown in Figure 16. The last digit 
may not be significant. 



model 


Mt 


a sw 






c 


K 


Gauss 


1/16 


0 


6.30 


6.30 


0.91 


0.333 




1/4 


0.01 


6.43 


6.44 


0.86 


0.34 




1 


0.25 


8.00 


8.40 


0.56 


0.46 




4 


0.66 


47.3 


86 


0.30 


0.50 


exponential 


1/16 


0 


6.38 


6.38 


0.62 


0.35 




1/4 


0 


6.89 


6.89 


0.68 


0.44 




1 


0 


13.85 


13.85 


0.46 


0.48 




4 


0 


339 


339 


0.28 


0.49 



As in the case of adhesive interactions, the normal displace- 
ment seems less sensitive to both the choice of the potential and 
the Tabor coefficient than the contact area, unless normal loads 
are very small, i.e., at loads similar in magnitude or smaller than 
the squeeze-out load for \ij = I. This is demonstrated in 
Figure 18. It reveals that information on the (effective) near- 
range surface interactions at small separation are difficult to 
obtain from experimentally measured load-displacement curves. 

I conclude this section with an analysis of the gap profile for 
repulsive contacts. At large loads, different Tabor parameters 
and functional forms for finite-range repulsion yield gap 
profiles that are indistinguishable at small magnification, see 
Figure 19a. Differences become nevertheless significant at high 
resolution near the center of the contact. Particularly remark- 
able is the data set for the Gauss model with [ij = 4 and its 
bistability revealed in Figure 19b. For an increasing force, no 
contact has formed at = 7.5. However, when reaching 
Fn = 7.5 from above, contact is formed for radii r < a c ~ 1.73. 
In the latter case, the gap then quickly increases within Ar ~ 0. 1 
to an almost constant value of order for r > a c , as if one 
had a single confined layer of liquid. For radii r > a™ cro , the 
gap assumes the "macroscopic" behavior. Here, a™ cm ~ 4 is 
the contact radius that one would ascertain from the analysis of 
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Figure 18: Displacement d as a function of normal force F N ' n the 
vicinity of the spontaneous wetting force F sw . Symbols reflect numeri- 
cal results. The lines, which connect many data points not explicitly 
shown, are drawn to guide the eye. The two thick grey lines reflect the 
square of the contact radius in the Hertz and DMT approximation, re- 
spectively. Color coding: uj = 4 (red), uj = 1 (green), and uj = 1/4 
(blue). 



the gap with low resolution, e.g., via graphical inspection of 
Figure 19a. 

At small loads, the sensitivity of the gap profile on the details of 
the model become even more apparent. This result, which can 
be seen in Figure 19c, is expected, since the elasticity of the tip 
is no longer relevant. Instead, the force-displacement curve is 
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Figure 19: Gap g(r) as a function of the lateral distance from the origin 
rfor a large load Fn = 60 (a) and (b) as well as for an intermediate load 
F N = 7.5 (c). In each case, surfaces repel each other. Graph (b) 
contains the same data as (a) but has higher resolution. Color coding: 
ut = 4 (red) and uj = 1 (green). 
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predominantly determined by the effective surface interactions, 
as shown clearly by the [ij = 4 data sets in Figure 18. They ex- 
hibit, to leading order, a oc exp(-J/Q relation in the expo- 
nential model, and a F N oc exp(-J 2 /^ 2 ) relation for the Gauss 
model, where £ is inversely proportional to 

Conclusion 

The principle new aspect of this work is the continuum- 
mechanics based analysis of single-asperity contacts with finite- 
range repulsion acting in addition to short-range hard-wall 
repulsion. The analysis is based on the concept of the Tabor 
coefficient and the repulsion is assumed to arise due to the pres- 
ence of a strongly wetting fluid. As for attractive single-asperity 
contacts, it is found that the contact area or the displacement on 
the normal load depend, to a large degree, not only on the 
surface energy but also on the Tabor coefficient jij- Moreover, 
for |i T exceeding a critical value, there may exist a range of 
loads in which two (meta)stable solutions coexist, i.e., one in 
which the surfaces touch and one in which a thin gap between 
the two surfaces remains. When the value for the load is 
increased above a threshold, the latter solution becomes 
unstable and the gap disappears. However, in order to obtain 
this kind of behavior, which is reminiscent of the squeeze-out of 
a wetting fluid, the finite-range interactions between the 
contacting surfaces have to be tailored correctly. Using a 
surface interaction Vf r , whose derivative increases monotoni- 
cally as the gap g approaches zero, such as Vf r oc exp(-g/zo), 
only one stable solution exists for any given normal load. 
Conversely, when the distance-force dependence is multi- 
valued, as is the case for a Vf r oc Qxp(-g 2 /2zq ^ relation, 
squeeze-out and spontaneous wetting can be rationalized and 
thus be modeled in the realm of continuum mechanics - in 
terms of transitions between (meta)stable solutions. These tran- 
sitions (similar to instabilities in the Prandtl model [26], in 
which a particle is dragged with a weak spring through a 
sinusoidal potential) can occur for solvated tips on surfaces, 
for example, if the effective tip-surface interactions has 
zero slope when the surfaces touch, as is the case for 



Vf r oc Qxp^-g 2 /2zq ^- In reality, the far-field potential may 
even be oscillatory and evidenced by the squeeze-out of many 
subsequent layers. Such behavior has been recently observed 
and linked to the (damped) long-range oscillatory behavior of 
the density correlations in high-density liquids [15,19]. 

An interesting consequence of short-range repulsion is that the 
contact geometry can look similar to that of an adhesive neck. 
This is shown in Figure 19b for the (jij = 4) Gauss model and 
decreasing load. To improve the visualization, a similar gap 
geometry is shown again in Figure 20 together with a profile of 
the finite-range repulsion. 

A secondary aspect of this work is devoted to the analysis of 
how to best reach well-defined asymptotic behavior in numeri- 
cal simulations of adhesive contact mechanics. It is found that 
the DMT limit is approached quickest when using attractive 
potentials whose first derivative disappears as the gap goes to 
zero, at least if the contact area is the variable of interest. 
However, these potentials approach the JKR limit only at a rate 
of l/jij for large |ij and the contact area becomes difficult to 
define once [ij>l. Thus, one is better off using potentials with 
finite slope in the small-gap limit. They converge in a well- 
defined fashion with \j juj to the JKR limit for large Tabor 
coefficients. This is supposedly the more relevant limit for 
adhesive surfaces with self-affine fractal roughness. For the 
modeling of repulsive surfaces, the situation is more compli- 
cated. Formally, the JKR limit is again reached more quickly 
with models that have finite slope at zero gap. However, these 
models do not allow one to model the hysteretic response of a 
confined fluid that results whenever the squeeze-out force 
exceeds the spontaneous wetting force. 

A by-product of this work is a minor modification of the 
phenomenological description of single-asperity contact 
mechanics by Carpick, Ogletree, and Salmeron [1]. The COS 
equations can be parametrized to contain the correct asymptotic 
behavior for JKR and for DMT limits and also for the superpo- 



^TTt^ ttTTV^ 



Figure 20: Contact geometry for a Gauss model with finite-range repulsion. Arrows indicate the direction of normal load Fn (thick arrow) and that of 
the finite-range repulsion (thin arrows) acting in addition to a hard-wall constraint. No adhesive forces between the surfaces are considered. 
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sition of extremely short- and long-range interfacial interac- 
tions, as shown by Schwarz [2]. However, they still have a few 
formal shortcomings for intermediate-range potentials. For 
example, the original interpolation of the contact-area-on-load 
dependence for finite Tabor coefficients recuperates neither 
Hertzian contact mechanics at large loads with correct prefac- 
tors nor the correct contact radii in either DMT or JKR limit at 
zero normal loads. In this work, I propose to enforce those 
limits exactly including the correct asymptotics for 
a c (\ij,F^ = 0) at |ix = 0 and jij — > 00 • By doing so, the 
maximum error of the a c (F N = 0,ur) curve could be reduced 
from 1.2% to less than 0.3%. A shortcoming of both the origi- 
nal and the new, modified COS equations is that they both 
assume an asymptotic behavior near pull-off (F^ — > -F p ) 
according to (a - a p ) oc (F p + F N ) K , where the exponent takes 
the JKR value k = 2/3 for any non-zero Tabor coefficient. The 
modified COS equations could thus be improved further if one 
incorporated the new finding that k crosses over continuously 
from 2/3 (exact for JKR) to 1/2 (exact for DMT). However, this 
does not seem useful in practice. Extreme accuracies (5 digits 
and more for a c and F N ) would be needed in measurements to 
deduce a p and k to within one or two digits. Such an accuracy is 
difficult to achieve both experimentally and numerically. More- 
over, the surface energy is not very well defined at small scales, 
because its precise value depends crucially on roughness down 
to the atomic scale, see, e.g., [27]. Thus, from a practical point 
of view, both the original and the modified COS equations are 
quite reasonable, all the more because the geometry of real tips 
can deviate quite substantially from a parabola. 

This work is concluded with an assessment of what values for 
jix one might expect in AFM or SFA experiments. To come up 
with a ballpark estimate, the following "typical values" shall be 
assumed: Ay = 40 mN (Ay can, of course, be close to zero, but 
much higher, e.g., for two equally charged surfaces in the 
context of electrochemistry), E = 5 GPa (in between soft matter 
and ceramics), zq = 10 A (size of an OMCTS or molten salt 
molecule), R = 1 um (in between AFM and SFA, precise value 
not very important, as third root is taken). These numbers lead 
to |ij = 0.4, which is close to the interesting "cross-over" 
regime. Thus, real contacts may span a broad range of values 
for jij- Comparison between theory and experiment may be 
difficult, in particular because atomic-scale roughness (or even 
sub-atomic roughness arising from electron orbitals) leads to 
complicated slip-boundary conditions and slow kinetics. 
However, given a well-motivated form for the effective inter- 
action between two flat surfaces, it may yet be possible to ratio- 
nalize and to model, at least on a semi-quantitative level, the 
interactions of curved surfaces in the presence of a strongly 
wetting fluid within the presented Tabor-coefficient based 
framework. Particularly appealing systems may be found in 



tribo-electrochemical applications, where the surface interac- 
tions can be tailored in a quasi-continuous fashion. 
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